use gdal::Dataset;
use gdal_sys::{GDALWarpAppOptionsNew, GDALWarp, GDALDatasetH};
use gdal::raster::RasterBand;
use image::{ImageBuffer, Rgba};
use std::fs;
use std::path::Path;
use std::error::Error;
use std::time::Instant;
use std::ffi::CString;
use std::ptr::null_mut;
use std::os::raw::c_char;

const TILE_SIZE: usize = 256;
const EARTH_CIRCUMFERENCE: f64 = 40075016.686;

fn warp_dataset(input_path: &str, output_path: &str) -> Result<Dataset, Box<dyn Error>> {
   
    let dataset = Dataset::open(input_path)?;

   
    let mut options = vec![
        CString::new("-of")?,
        CString::new("MEM")?,
        CString::new("-r")?,
        CString::new("near")?,  
        CString::new("-t_srs")?,
        CString::new("EPSG:3857")?,
        CString::new("-wo")?,
        CString::new("NUM_THREADS=40")?,
    ];

   
    let mut c_options: Vec<*mut c_char> = options
        .iter_mut()
        .map(|s| s.as_ptr() as *mut c_char)
        .collect();


  
    let dest_path = CString::new(output_path)?;

    unsafe {
        let warp_options = GDALWarpAppOptionsNew(c_options.as_mut_ptr(), null_mut());
        if warp_options.is_null() {
            return Err("Failed to create warp options".into());
        }

        let mut c_dataset: GDALDatasetH = dataset.c_dataset();
        let dst_dataset = GDALWarp(
            dest_path.as_ptr(),
            null_mut(),
            1,
            &mut c_dataset,
            warp_options,
            null_mut()
        );

        if dst_dataset.is_null() {
            return Err("Warp operation failed".into());
        }

        let _res_dataset = Dataset::from_c_dataset(dst_dataset);

        Ok(_res_dataset)
    }
}

fn web_mercator_to_tile(x: f64, y: f64, zoom: u32) -> (u32, u32) {
    let n = 2.0f64.powi(zoom as i32);
    let xtile = ((x + EARTH_CIRCUMFERENCE / 2.0) * n / EARTH_CIRCUMFERENCE).floor() as u32;
    let ytile = ((EARTH_CIRCUMFERENCE / 2.0 - y) * n / EARTH_CIRCUMFERENCE).floor() as u32;
    (xtile, ytile)
}

fn tile_to_web_mercator(xtile: u32, ytile: u32, zoom: u32) -> (f64, f64) {
    let n = 2.0f64.powi(zoom as i32);
    let x = xtile as f64 * EARTH_CIRCUMFERENCE / n - EARTH_CIRCUMFERENCE / 2.0;
    let y = EARTH_CIRCUMFERENCE / 2.0 - ytile as f64 * EARTH_CIRCUMFERENCE / n;
    (x, y)
}

fn generate_tile(dataset: &Dataset, zoom: u32, tx: u32, ty: u32, tile_path: &Path) -> Result<(), Box<dyn Error>> {
    let (tile_x_min, tile_y_max) = tile_to_web_mercator(tx, ty, zoom);
    let (tile_x_max, tile_y_min) = tile_to_web_mercator(tx + 1, ty + 1, zoom);

    let mut img = ImageBuffer::new(TILE_SIZE as u32, TILE_SIZE as u32);

    let transform = dataset.geo_transform()?;
    let (origin_x, pixel_width, _, origin_y, _, pixel_height) = (
        transform[0], transform[1], transform[2],
        transform[3], transform[4], transform[5]
    );

    let raster_size = dataset.raster_size();

    let band1: RasterBand = dataset.rasterband(1)?;
    let band2: RasterBand = dataset.rasterband(2)?;
    let band3: RasterBand = dataset.rasterband(3)?;
    
    let nodata1 = band1.no_data_value().unwrap_or(f64::NAN);
    let nodata2 = band2.no_data_value().unwrap_or(f64::NAN);
    let nodata3 = band3.no_data_value().unwrap_or(f64::NAN);

    for (x, y, pixel) in img.enumerate_pixels_mut() {
        let merc_x = tile_x_min + (tile_x_max - tile_x_min) * (x as f64 / TILE_SIZE as f64);
        let merc_y = tile_y_max - (tile_y_max - tile_y_min) * (y as f64 / TILE_SIZE as f64);

        let px = (merc_x - origin_x) / pixel_width;
        let py = (merc_y - origin_y) / pixel_height;

        if px >= 0.0 && px < raster_size.0 as f64 && py >= 0.0 && py < raster_size.1 as f64 {
            let (red, green, blue) = if zoom <= 3 {
                let (avg_red, avg_green, avg_blue) = compute_average_color(&band1, &band2, &band3, px, py, raster_size, 24, nodata1, nodata2, nodata3)?;
                (avg_red, avg_green, avg_blue)
            } 
            else if zoom <= 5 {
                let (avg_red, avg_green, avg_blue) = compute_average_color(&band1, &band2, &band3, px, py, raster_size, 18, nodata1, nodata2, nodata3)?;
                (avg_red, avg_green, avg_blue)
            } 

            else if zoom <= 7 {
                let (avg_red, avg_green, avg_blue) = compute_average_color(&band1, &band2, &band3, px, py, raster_size, 16, nodata1, nodata2, nodata3)?;
                (avg_red, avg_green, avg_blue)
            } 

            else if zoom <= 11 {
                let (avg_red, avg_green, avg_blue) = compute_average_color(&band1, &band2, &band3, px, py, raster_size, 10, nodata1, nodata2, nodata3)?;
                (avg_red, avg_green, avg_blue)
            }

            else if zoom <= 15 {
                let (avg_red, avg_green, avg_blue) = compute_average_color(&band1, &band2, &band3, px, py, raster_size, 6, nodata1, nodata2, nodata3)?;
                (avg_red, avg_green, avg_blue)
            }
            else if zoom <= 17 {
                let (avg_red, avg_green, avg_blue) = compute_average_color(&band1, &band2, &band3, px, py, raster_size, 2, nodata1, nodata2, nodata3)?;
                (avg_red, avg_green, avg_blue)
            } 
            else {
                let red = band1.read_as::<f64>((px as isize, py as isize), (1, 1), (1, 1), None)?;
                let green = band2.read_as::<f64>((px as isize, py as isize), (1, 1), (1, 1), None)?;
                let blue = band3.read_as::<f64>((px as isize, py as isize), (1, 1), (1, 1), None)?;
                (red.data[0], green.data[0], blue.data[0])
            };

            if red == nodata1 || green == nodata2 || blue == nodata3 {
                *pixel = Rgba([0, 0, 0, 0]); 
            } else {
                *pixel = Rgba([
                    red.clamp(0.0, 255.0) as u8,
                    green.clamp(0.0, 255.0) as u8,
                    blue.clamp(0.0, 255.0) as u8,
                    255
                ]);
            }
        } else {
            *pixel = Rgba([0, 0, 0, 0]);
        }
    }

    img.save(tile_path)?;

    Ok(())
}

fn compute_average_color( band1: &RasterBand, band2: &RasterBand, band3: &RasterBand,px: f64, py: f64, raster_size: (usize, usize),window_size: isize, nodata1: f64, nodata2: f64, nodata3: f64) -> Result<(f64, f64, f64), Box<dyn Error>> 
{
    let px_int = px.floor() as isize;
    let py_int = py.floor() as isize;
    let x_offset = px_int.max(0).min(raster_size.0 as isize - window_size);
    let y_offset = py_int.max(0).min(raster_size.1 as isize - window_size);
    let width = (raster_size.0 as isize - x_offset).min(window_size) as usize;
    let height = (raster_size.1 as isize - y_offset).min(window_size) as usize;

    let red = band1.read_as::<f64>((x_offset, y_offset), (width, height), (width, height), None)?;
    let green = band2.read_as::<f64>((x_offset, y_offset), (width, height), (width, height), None)?;
    let blue = band3.read_as::<f64>((x_offset, y_offset), (width, height), (width, height), None)?;

    if red.data.iter().any(|&v| v == nodata1) || green.data.iter().any(|&v| v == nodata2) || blue.data.iter().any(|&v| v == nodata3) {
        return Ok((nodata1, nodata2, nodata3));
    }

    let pixel_count = width * height;
    let avg_red = red.data.iter().sum::<f64>() / pixel_count as f64;
    let avg_green = green.data.iter().sum::<f64>() / pixel_count as f64;
    let avg_blue = blue.data.iter().sum::<f64>() / pixel_count as f64;

    Ok((avg_red, avg_green, avg_blue))
}


fn main() -> Result<(), Box<dyn Error>> {
    let start_time = Instant::now();

   
    let input_path=r"C:\Users\DELL\Desktop\RD0100_DataCloud_M52D008008_2024H1REP_4FF.tif";
    let output_path=r"C:\Users\DELL\Desktop\bbb_3857.tif";
   
    let warped_dataset = warp_dataset(&input_path, output_path)?;

    let geo_transform = warped_dataset.geo_transform()?;
    let (origin_x, pixel_width, _, origin_y, _, pixel_height) = (
        geo_transform[0], geo_transform[1], geo_transform[2],
        geo_transform[3], geo_transform[4], geo_transform[5]
    );

    let raster_size = warped_dataset.raster_size();
    let (raster_width, raster_height) = (raster_size.0 as f64, raster_size.1 as f64);

    let x_min = origin_x;
    let y_min = origin_y + pixel_height * raster_height;
    let x_max = origin_x + pixel_width * raster_width;
    let y_max = origin_y;

    for zoom in (0..=17).rev() {
        let (tx_min, ty_max) = web_mercator_to_tile(x_min, y_min, zoom);
        let (tx_max, ty_min) = web_mercator_to_tile(x_max, y_max, zoom);

        for ty in ty_min..=ty_max {
            for tx in tx_min..=tx_max {
                let tile_path_str = format!(r"C:\Users\DELL\Desktop\test/{}/{}/{}.png", zoom, tx, ty);
                let tile_path = Path::new(&tile_path_str);
                fs::create_dir_all(tile_path.parent().unwrap())?;
                generate_tile(&warped_dataset, zoom, tx, ty, tile_path)?;
            }
        }
    }
    let duration = start_time.elapsed();
    println!("TotalTime: {:?}", duration);
    Ok(())
}
